PASJ: Publ. Astron. Soc. Japan 59, 1-??, 2007 October 25 
© 2007. Astronomical Society of Japan. 



Self-gravitational Magnetohydrodynamics with Adaptive Mesh 
Refinement for Protostellar Collapse 

Tomoaki Matsumoto 

Faculty of Humanity and Environment, Hosei University, Fujimi, Chiyoda-ku, Tokyo 102-8160, Japan 

matsu@i. hosei.ac.jp 

(Received 2006 September 5; accepted 2007 June 25) 
Abstract 

A new numerical code, called SFUMATO, for solving self-gravitational magnetohydrodynamics (MHD) 
problems using adaptive mesh refinement (AMR) is presented. A block-structured grid is adopted as the 
grid of the AMR hierarchy. The total variation diminishing (TVD) cell-centered scheme is adopted as 
the MHD solver, with hyperbolic cleaning of the divergence error of the magnetic field also implemented. 
The MHD solver exhibits a second-order accuracy in convergence tests of linearized MHD waves. The 
self-gravity is solved using a multigrid method composed of (1) a full multigrid (FMG)-cycle on the AMR 
hierarchical grids. (2) a V-cycle on these grids, and (3) an FMG-cycle on the base grid. The multigrid 
method exhibits spatial second-order accuracy, fast convergence, and scalability. The numerical fluxes are 
conserved by using a refluxing procedure in both the MHD solver and the multigrid method. Several tests 
are performed and the results indicate that the solutions are consistent with previously published results. 
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1. Introduction 

Protostellar collapse is one of the most important pro- 
cesses in star formation; it exhibits a high dynamic range 
in density and spatial dimensions. Adaptive mesh refine- 
ment (AMR) is a powerful technique for performing nu- 
merical simulations requiring a high dynamic range in the 
mesh schemes. Local high-resolution is realized by em- 
ploying grids of differing resolutions. The finer grids are 
inserted and their location changed according to given re- 
finement criteria. Following the introduction of the fun- 
damental principles of AMR by Berger & Oliger (1984) 
and Berger & Colclla (1989), this technique is becoming 
widely used in astrophysical simulations. 

In star formation, along with self-gravity, the mag- 
netic field also plays an important role. Both self-gravity 
and magnetohydrodynamics (MHD) are therefore imple- 
mented in recent AMR codes (see the review by Klein et 
al. 2006). 

Existing self-gravitational MHD AMR codes are gener- 
ally based on Cartesian grids; although block-structured 
grids tend to be most often frequently, alternative ap- 
proaches also exist, e.g. RAMSES (Fromang et al. 2006). 
In block-structured AMR, numerical cells arc refined in 
the unit of the block, and a block is itself an ordinal uni- 
form grid. Block-structured grids are divided into two 
categories: patch-oriented grids and self-similar blocks. 
In the former approach, a block has variable number of 
cells, and the shape of the block is also changed in the 
course of refinement. A large block containing many cells 
and a small block containing just a few cells can co-exist. 
This approach originated in the study of Berger & Colella 
(1989), and has been adopted in many codes: e.g., the 



AMR code Orion of Berkeley (Truelove et al. 1998), Enzo 
(Norman & Bryan 1999), and RIEMANN (Balsara 2001). 
In the latter approach, all blocks contain the same num- 
ber of cells, but the physical sizes of the blocks are differ- 
ent. This approach has also been adopted by many codes: 
e.g., FLASH (Fryxell et al. 2000; Banerjee & Pudritz 
2006), and NIRVANA (Ziegler 2005). This approach of- 
fers advantages including relatively simple algorithms for 
refinement, parallelization, and vectorization. This paper 
also follows this latter approach. In particular, this ap- 
proach enables the AMR code to implement a full multi- 
grid (FMG)-cycle in the multigrid method for self-gravity. 
This scheme is an extension of the multigrid method for 
the nested grid approach (Matsumoto & Hanawa 2003a) 
to the AMR grid. 

There are several ways to treat the MHD, particularly in 
the treatment of the V • B term. These approaches may 
be categorized as (1) the constrained transport method 
with a staggered grid, (2) the projection scheme with a 
Poisson solver, and (3) the eight-wave formulation (see 
Toth 2000; Balsara & Kim 2004 for a comparison of these 
schemes) . 

Recently, Dedner et al. (2002) proposed an alternative 
scheme for cleaning the V ■ B term. This scheme has a 
formulation that is similar to, but distinct from, the eight- 
wave formulation; two additional waves transfer the V • B 
error isotropically at a given speed, and the waves decay at 
a given rate. As a result, the error is propagated indepen- 
dently of the gas motion, and so is diluted. By contrast, 
using the eight- wave formulation, a single additional wave 
transfers the error at the gas velocity. The error is always 
propagated downstream of the gas motion, and can be- 
come stagnated at the shock wave and the center of the 
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collapsing cloud. The scheme of Dedner et al. (2002) is 
therefore adopted here. 

The author has previously developed a self-gravitational 
MHD code with a nested grid (Matsumoto & Hanawa 
2003a; Matsumoto & Hanawa 2003b; Matsumoto et al. 
2004; Machida et al. 2004). In the nested grid approach, 
the grids are refined adaptively, while the sizes and posi- 
tions of the sub-grids are fixed. The problems to which 
this method may be applied have therefore been restricted 
to those problems whose objects that require fine grid res- 
olution are located in the center of the computational do- 
main, e.g., Machida et al. (2006). 

A self-gravitational MHD code with AMR has now, 
however, been developed, and this method is presented 
in this paper. In § 2-4, the governing equations to be 
solved, the rules of discretization, and the refinement al- 
gorithm are presented. In § 5, the present methods for 
treating the hydrodynamics and MHD of the problem are 
described, while in § 6, the methods of the multigrid it- 
eration for solving self-gravity are described. In § 8, the 
results of several numerical tests are presented. The paper 
is concluded in § 9. 

2. Governing Equations 

The equations of self-gravitational, ideal hydrodynam- 
ics and MHD can be expressed in a conservative form with 
a source term, 



dU 
~dt 



dF x dF y dF z 



= S, 



(1) 



dx dy dz 

together with Poisson's equation 

V 2 $ = AirGp, (2) 

where U is a vector of conservative variables, F x , F y , and 
F z are numerical fluxes, and S is the source term vec- 
tor. In equation (2), <f>, G, and p denote the gravitational 
potential, the gravitational constant, and the density, re- 
spectively. 

For ideal MHD with self-gravity, the vectors in equa- 
tion (1) are expressed by 



U = (p,pv x ,pv y ,pv z ,B x ,B y ,B z ,pE) 



( 



P v l 



f P+i 

pV x Vy 

pv x v z 



pv x 

B\ 2 /St: -B 2 x /4n 



\ 



(3) 



(4) 



B x B y /4ir 
B x B z /4n 


^xBy VyB x 

~v z B x +v x B z 
V (pE + P+\B\ 2 /8ir)v x -B x {B-v)/An J 



S= (0,P9x,P9 y ,P9z, 0,0,0, pg-v) , (5) 

where v = (v x , v y , v z ) T represents velocity, B = 
(B x , By, B Z ) T represents the magnetic field, g = 
(9x,g y ,9z) T = -V$ represents gravity, E = \ v\ 2 /2 + (7 - 
l)~ 1 P/p + \B\ 2 /8iTp is the total energy, and P represents 



pressure. The vectors F y and F z are obtained by rotating 
the components in F x by the right-hand rule. 

For ideal hydrodynamics with self-gravity, the govern- 
ing equations arc obtained by setting B — and omitting 
the 5th to the 7th components of equation (1) and vectors 
(3)-(5). 

Barotropic and isothermal equations of state are also 
implemented in the AMR code. In these equations of 
state, the component corresponding to E in equation (1) 
and vectors (3)-(5) (the eighth component) are excluded, 
and P is expressed as a function of p. 

3. Discretization 

Equations (1) and (2) are solved by a difference scheme 
based on the finite-volume approach. The computational 
domain is divided into cells, each of size Ax x Ay x Az. 
Each cell is labeled by (i,j,k), the indices of the cell 
in the x, y, and z-directions, respectively. The loca- 
tion of the cell center is indicated by the position vec- 
tor r.ij.k- The conservative variables U, the source term 
S, and the gravitational potential $ are defined at the 
cell center, i.e., U itj . k :=U{r it j^), S itj , k := S(r it j t k), and 



4> 



i,3,k 



j,k)- 



The numerical fluxes F x , F v , and F z 



arc defined at the cell surfaces with normals in the x, y, 
and z-directions, respectively. For convenience, the nota- 
tion of F x ^ i±1 / 2 j.k '■— Fx [ r ij,k ± (Ax/2)x] is introduced, 
where x denotes the unit vector in the x-direction. We 
introduce the following notation to describe the spatial 
differences, 



9xQi+l/2,j,k — 
3 X Qi.j : k 



Qi+lJ.k Qi,j,k 

Ax ' 

Qi+l/2,j,k — Qi-l/2,j,k 

Ax ' 



q2 n _ Qi+l.j,k — 2Qi,j.k + Qi-l,j,k 



(6) 
(7) 
(8) 



The differences in the y- and z-dircctions arc expressed in 
a similar manner. 

4. Grid Refinement 

A self-similar block-structured grid is adopted. Each 
block consists of N x x N y x N z cells, where N x , N y , and 
N z denote the number of cells in the x, y, and z-directions 
respectively, with all blocks having the same number of 
cells. The number of cells inside a block is fixed in course 
of calculation, but the cell width differs depending on the 
grid-level. The blocks are thus self-similar. A schematic 
diagram of a block-structured grid is shown in Figure 1. 
The number of cells is set at N x = N y = N z = 8 in this 
figure. 

If some cells satisfy a refinement criterion, the block in 
which these cells lie is divided into 8 child blocks, and 
every cell inside the parent block is also refined into 8 
child cells (see Fig. 1). The cell width of the child cells is 
half of the parent cell width, that is the refinement ratio is 
fixed at two in this implementation (in some AMRs having 
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Fig. 1. Schematic diagram of the grid refinement process. 
The thick lines represent the boundaries of the blocks and the 
thin lines represent the boundaries of the cells. Each block 
consists of 8 3 cells in this figure. The hatched block is refined 
into 8 child blocks. 

patch-oriented grids any refinement ratio of the power of 2 
can be used). The coarsest grid-level is labeled £ — (the 
base grid), and the finest grid- level is labeled I — £ max . 
The £-th grid-level has a 2 times higher spatial resolution 
than the coarsest grid-level. The block-structured grid is 
managed by an octree structure; the parent (coarse) block 
is linked with eight fine (child) blocks. Moreover, a block 
is linked with its neighboring blocks. These link lists are 
reconstructed every time the grids are refined. 

In the construction of a child block, the conservative 
variables of the parent cells Uf JK are interpolated to 



obtain those of the child cells U 



For this purpose 



linear interpolation with a slope limiter is adopted, 



Utt , = Uf lir + VU" ■ (rt, fc - rfr K ), 



(9) 



where r^j k 



i,j,k — ^ I,J,K t v ^ ' V i,j,k 

and rf j K denote the position vectors indicat- 
ing the centers of the child and parent cells, respectively. 
The gradient inside the parent cell is slope-limited accord- 
ing to 



/ minmod (d* £7^,2 ,j !K ,d x Uf_ 1/2 JK ) \ 



vu H = 



minmod [dyUf >J+1/2tK ,dyU It j_ 1/2iK 
y minmod (d z Uf JK+1/2 ,d z U? ilK _ 1/2 



2. If blocks of grid- level 1 + 2 exist, then the corre- 
sponding blocks of t are also marked. 

3. Blocks adjacent to marked blocks are also marked, 
so that the grids of level £+1 properly nest the grids 
of level £ + 2. 

4. New blocks of level £ + 1 are constructed as child 
blocks of the marked blocks. 

5. Blocks of £+ 1 are removed if their parent blocks are 
not marked. 

These procedures are called in ascending order of grid- 
level, from £ max — 1 to 0. 

5. Hydrodynamics and MHD 

5.1. Basic Solvers 

The governing MHD equations (1), which include as 
a special case the hydrodynamics equations, are solved 
on the block-structured grid described in § 4. Methods 
used on ordinal uniform grids can be applied to the 
block-structured grid if boundary conditions are properly 
specified for each block. A monotone upstream-centered 
scheme for conservation laws (MUSCL) approach and 
predictor-corrector method are adopted here for integra- 
tion with respect to time in order to achieve a second- 
order accuracy in space and time (e.g., Hirsch 1990), and 
an unsplit approach rather than a a fractional timcstep 
approach is adopted. 

The numerical flux is obtained using the linearized 
Riemann solver; the solvers are based on the schemes of 
Roe (1981) and Fukuda & Hanawa (1999). For MHD, the 
hyperbolic divergence cleaning of Dedncr ct al. (2002) is 
adopted for reducing V • B. According to Dedner et al. 
(2002), equations (3)-(5) are then modified as follows: 

U = (p 7 pv X7 pv y ,pv z ,B X7 B y ,B Z7 pE,ifj) T , (12) 



,(10) 



/ 



where minmod(-,-) denotes the minmod function. Note 
that this interpolation conserves the conservative variables 
in the refinement procedure, 



/ PVx 

pvl + P+\B\ 2 /8n-Bl/4n 

pV X Vy - BxBy/4-K 

pv x v z -B x B z /Att 

V X By - VyB X 

-v z B x + v x B z 
(pE + P+ \B\ 2 /8tt)v x - B X (B ■ v)/4tt 
clB x 



\ 



,(13) 



la' 



U H {r)dr 



U h {r)dr, 



(11) 



where Qfj K denotes the zone of a parent cell whose cen- 
ter is located at rfj K . 

The refinement algorithm is based on that of Berger & 
Colella (1989), where grids of level £ are refined using the 
following procedures to construct grids of level £+1: 

1. A block of grid-level £ is marked if the cells inside 
the block satisfy a refinement criterion. 







P9x 
P9y 
P9z 



(V-B)B X /4:TT 
{V-B)B y /A-K 
{\7-B)B z /4n 

(14) 



•«)-B-(V^)/4tt 

V -(c h /c p ) 2 ip 

where ip is a scalar potential propagating a divergence 
error, Ch is the wave speed, and c p is the damping rate 
of the wave. Note that this modification requires addi- 
tional components in the basic equation, and a total of 
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nine waves are solved for. This approach is similar to 
the eight-wave formulation. In the method of Dcdner et 
al. (2002), two additional waves transfer the V • B error 
isotropically at a speed of c/j, and the waves decay at a 
rate of c p . Similarly, in the eight-wave formulation, an 
additional wave transfers the divergence error at a gas 
velocity of v. 

In the predictor step, U™^ k is updated to U™^ 2 
half time step, i.e., 



by a 



U 



n+l/2 
i,3,k 



u i.j,k 



x,i,j,k 



(15) 



where the superscript n denotes the time level, and At = 
t n+1 — t n . The numerical flux F n , which is defined at the 
cell boundaries, is calculated using the primitive variables, 

Q= {p,v x ,v y ,v z ,B x ,By,B z ,P,il>) T . (16) 

The numerical flux has second-order spatial accuracy due 
to the MUSCL extrapolation of the amplitudes of the 
eigenmode of the MHD waves, LU, where L denotes the 
matrix of the left eigenvector (see the appendix in Fukuda 
& Hanawa 1999 for details) 1 . 



The source term S i \ k is calculated using 



( 



n n-1/2 
Pi,j,k9 x .i,j,k 

n n-1/2 
Pi,j,k9y^i^j^k 

n n-1/2 
Pi,j,k9 z,i,j,k 



(17) 





(V-B 

(V-B) 

(V-B 









/4tt 



n jjn 

i,j,k x,i,j,k 



\ 



n C "-V2 

Pi,j,k\ y i,j,k 







(V^)™,, fc /47T 



/ 



where (V • B) itj , k = d x B x ^^_ k + d y B y ^j^ k + d z B z ^j^ k , and 
the magnetic field defined at the cell surfaces, B x i ±i/ 2 j. k , 
is obtained from the ninth component of F x ,i±i/2,j,k- 

= (d x ipi,j,k, dytpij,!,, dz^ij,k) T is 
fifth to seventh components of 
The ninth com- 
ponent of the source term (14) is evaluated separately by 
operator splitting, in which the formal solution is used as 
follows, 



Similarly, (VV>)ij,fc 
obtained from the 

F x,i±l/2,j,ki F y,i,j±l/2,ki an d F z ,iJ y k±l/2 







t+1/2 _ 



-1/2, 



exp 



At 
~2~ 



(18) 



where ^™+ 1 / 2 ,™* i s the ninth component of U n+1 / 2 solved 
by equation (15). The free parameters cu and c p are re- 
lated to the time-marching, and are described in § 5.2. 



1 For the case of the hydrodynamics, the MUSCL extrapolation 
is applied to the primitive variables, Q in the predictor and 
corrector steps in this implementation. 



Note that gravity g™~ k lags by half a time step. This 
slight lagging is expected to have a negligible effect on 
the accuracy in the predictor step, and the gravity in the 
previous corrector step can be reused in this predictor step 
(Truelove et al. 1998). This avoids an additional call of the 
multigrid method to solve Poisson's equation, significantly 
reducing the computational costs of this method. 

In the corrector step, a spatially second-order numeri- 
cal flux F n+1/2 is obtained by applying MUSCL extrap- 
olation to the amplitudes of the eigenmodes, which are 
converted from J7 n+1 / 2 . Using this flux, U^j k is updated 
to U™~j k by a full time step, 

U i,j,k — U i,j,k 



' '*■> F'j. ]j k 



g pn+1/2 
u V r y,i,j,k 



g F n+l/2 

" zf z,i,j,k 



? n+l/2\ 
} i,j,k J 



(19) 



The source term in the corrector step is estimated at the 
time leveH = i" +1 / 2 , 



qn+l/2 

i,j,k 



n+l/2 n+l/2 
Pi,j,k 9x.i.j,k 

n + l/2 n+l/2 
Pi,j,k 9y,i,j,k 

n+l/2 n+l/2 
PiJ,k @z,i,j,k 





(V • B)?+l /2 B 



x,i,j,k /4 71 " 

B:+V?/4n 



(20) 
\ 



'i,j,k 

V *'i,j-k ^y^ijjjk 

(V-B)^ 2 <tff/4' 






V 



n+l/2 
Pi,3,k 



(a n+1 J 2 



n+l/2 
'i,j,k 







R n+l/2 
°i,3,k 



where p™~^]/ 2 and v^]/ 2 are obtained from U™^ 2 , and 



'i,3,k 



i,3,k 



g"^ k is obtained by solving Poisson's equation (2) at the 
half time step as follows, 



i,j,k 



1/2 , n n+l/2 



(21) 



Poisson's equation is solved by the multigrid method, as 
shown in § 6, and the multigrid method prepares not only 
the cell-centered $, but also the cell-surfaced g. Gravity 
at the cell center g™^ 2 is obtained by averaging the val- 
ues of gravity at the cell surfaces, 

/ 

1 
2 



n+l/2 



n+l/2 n+l/2 
9 x ,i+l/2,j,k 9x,i-l/2,j,k 
1+1/2 n+l/2 



\ 



1y,i,j-l/2,k 



9y,i,j + l/2,k 
n+l/2 n+l/2 
\ 9z.i.j,k+l/2 + 9z,i.j,k-l/2 ) 



(22) 



Note that g at a cell surface is obtained by the dif- 
ference of the cell-centered $ (see equations [42]~[44]). 
Equation (22) therefore coincides with the central differ- 
ence of $ as far as the inside of a block is concerned. 

Similar to the predictor step, the ninth component of 
the source term is evaluated by operator splitting, 

2" 

/ C i, 

exp — At I ■ 



n+l 



n+l,* 



(23) 
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Fig. 2. Adaptive time-step scheme. The thick arrows de- 
note time steps, while the associated numbers, Q, ©, @, 
indicate the order of time-marching. 

where ^ n+1 >* is the ninth component of U n+1 obtained 
by equation (19). 

5.2. Time Marching 

The time-marching represented by equations (15) and 
(19) proceeds in units at the grid-level. The code is 
equipped with two modes of time-marching: an adaptive 
and a synchronous time-step mode. Which mode is used 
depends on whether the gas is non-self-gravitational or 
self-gravitational. The adaptive time-step mode is appro- 
priate for non-self-gravitational gases, and is based on the 
method of Berger & Colella (1989). In this mode, a coarser 
grid has a longer time step than a finer grid. In contrast, 
the synchronous mode is appropriate for self-gravitational 
gases, and every grid-level has the same time step. This 
is because evolution on the fine grid affects the detached 
coarse grid immediately due to self-gravity, and so the 
same time step must to be chosen for every grid-level (see 
discussion of Truelove et al. 1998). Note that the adaptive 
time-step mode could also be used for a self-gravitational 
gas at the expense of the first-order temporal accuracy of 
the scheme. 

Figure 2 shows the order in which the grid-levels pro- 
ceed for the adaptive time-step mode schematically. The 
numbers associated with the thick arrows denote the order 
of time-marching. Coarser grid-levels precede finer grid- 
levels. The fine grid-level undergoes several sub-cycles 
until the time level of the fine grid-level is synchronized 
with that of the coarse grid-level. The time step of a finer 
grid-level, At h , is given by 

At h = At H 2~ n , 

n = min {m £ N | At H 2~ m < A^ FL } . 

where Af^ FL denotes the time step calculated directly by 
the CFL condition at the fine grid-level, and At H denotes 
the time step at the coarser grid-level. Note that At h 
is fixed at At H /r in Berger & Colella (1989), where r 
denotes a refinement ratio. On the other hand, in the 
method presented here At can be equal to At H /2 n for 
n = 0, 1,2, if At' 1 satisfies the CFL condition. In the 
synchronous time-step mode, a common time step At is 



(24) 
(25) 



used at all grid- levels, and is given by 
At= min \At e rvj } , 



(26) 



where At^ FL denotes the time step calculated by the CFL 
condition at grid-level I. 

The solver described in § 5.1 includes two parameters, 
Ch and c p , and these are related to the mode of the time- 
marching. The wave speed Ch is obtained as 

ch = CFL J^, (27) 

where CFL denotes the CFL number, and At e=0 denotes 
the time step at grid-level £ = 0. For the adaptive time- 
step mode, 

(28) 



h = mm{Ax,Ay,Az e } 



while for the synchronous time-step mode, 
h= min {Ax e ,Ay e ,Az e } , 



(29) 



In both cases, Ch is constant across all the grid-levels, and 
satisfies the CFL condition at every grid-level. 
The damping rate c p is obtained from Ch, 



0.18Lc h , 



(30) 



where L is a scale length of a problem, and the coefficient 
of 0.18 is chosen according to Dcdner et al. (2002). This 
coefficient is a free parameter specifying the ratio of the 
damping rate and propagation speed of V • B. We con- 
firmed that this selection worked well in many calculations 
(e.g., Machida et al. 2005a; Machida et al. 2005b) and in 
the numerical tests presented in this paper. 

5.3. Boundary Condition for Ghost Cells 

Each block has ghost cells overlapping with the adjacent 
blocks. Figure 3 shows the interface between the coarse 
and fine blocks in the cc-direction. Two ghost cells are 
prepared in each direction, as indicated by the gray circles 
in the figure, since the MUSCL extrapolation requires two 
cells on both the left and right sides of a cell boundary 
where numerical flux is obtained. 

A boundary condition is imposed on the ghost cells 
each time before the predictor and corrector steps pro- 
ceed. When a block is adjacent to another block of the 
same grid-level, data is exchanged between the two blocks 
by a simple copy. When a block is adjacent to a coarse 
block, data is interpolated spatially. In addition, data is 
interpolated temporally if an adaptive time-step mode is 
used. The procedure for the boundary conditions is as 
follows: 

1. Adjacent coarse blocks are identified. 

2. Conservative variables on the coarse cells overlap- 
ping the ghost cells of a fine block arc interpolated 
temporally to coincide with the time level of the 
fine grid-level. This procedure is omitted if the syn- 
chronous time-step mode is used. 
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Fig. 3. The interface between a coarse block and a fine block. The grid points of the ghost cells for the fine block are denoted by 
gray circles, while the real grid points are denoted by filled circles. The arrows denote the numerical fluxes through the interface 
between the fine and coarse blocks. 



3. The temporally interpolated variables are converted 
into primitive variables. They are interpolated spa- 
tially to coincide with the grid points of the ghost 
cells of the fine grid-level. The interpolated variables 
are copied into the ghost cells. 

4. Adjacent blocks in the same grid-level arc identified. 

5. The variables of cells overlapping the ghost cells are 
simply copied into the ghost cells. 

In the temporal interpolation, the quadratic interpola- 



H,n 



U H ' n+1 '\ and U H > n+1 . Note 



tion is performed on U 
that U H ' n and U H > n+1 have of second-order temporal ac- 
curacy while {j H ' n+1 / 2 i s f first-order temporal accuracy. 
Using these variables, U H (t) in t n < t < t n+1 is interpo- 
lated as, 



U H {t) = (u 



T H,n+l 



t-t n 



t)U 



t n+1 - t 

+ (t-t n )U (lhH ' n+1 



(31) 

where lj( 1 )> H ' n+1 denotes the conservative variables of 
first-order accuracy at the time level t n+1 , defined as, 

jj{l),H,n+l _ 2jjH,n+l/2 _ jjH.n (32) 

A slope-limited gradient is adopted for the spatial in- 
terpolation. In Figure 3, the gradient of the primitive 
variables inside the cell (I,J,K) is evaluated using 



VQ 



H 



d x Qf-l/2,J,K 

minmod (dyQ"j +1/2 K ,d y Q 1 j I j_ 



minmod 



1/2, K 
l/2i®zQ?,J,K-l/2 



It may be noted that equation (33) describes a simple 
interpolation in the normal direction, with a slope-limited 
interpolation in the transverse directions. 

5.4- Refluxing at the Interfaces 

Rcfluxing was introduced by Berger & Colella (1989), 
and has often been used in standard implementations of 
AMRs. This technique maintains consistency between the 
fine and coarse grid-levels during time-marching, and en- 
sures that conservation laws, e.g., mass conservation, are 
satisfied. 

As shown in Figure 3, Uf JK is updated using 
F x,i+i/2,J,k, while U h t3 % is updated using F^"_ 1/2 j >s 
through several sub-cycles until the fine grid-level catches 
up with the coarse grid-level, where j G + 1], k(z[k,k + 
1], and n denotes the index of the sub-cycles of the fine 
grid-level. Obviously, F" I+1/2 j K AS H At H should equal 
£ F , , - AS h At h ' n to ensure that the conserva- 

j,K,n x,i— 1/2, j.k 

tion laws are satisfied, where AS H = Ay H Az H and AS*' 1 = 
Ay h Az h . In rcfluxing, Uf j K is re-calculated, taking ac- 



count of the difference between F 



j,k,ti x a. 



, . ~AS h At h ' n . 

-l/2,J,fe 



H 

x, 1+1/2, J,K 



AS H At H and 



More explicitly, as indicated 



in Figure 4, the refluxing procedure updates U 
according to 



H. 



to U 



H 



U H = U H <* 

1 



AV H 



.(33) 



E 

n surface 



F h x ' n At h ' n AS h 



F"At H AS 



H 



(34) 



where 
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Fig. 4. Rcfluxing procedure at the interface between fine 
and coarse cells. 

n 

J2 AS h = AS H , (36) 

surface 

and AV H denotes the volume of the coarse cell. 

6. Multigrid Method for Self-gravity 

6.1. Multigrid Cycles 

The multigrid method is widely used in many AMRs 
for solving Poisson's equation. In many AMRs, a solu- 
tion converges by means of the V-cycle on the grids of 
AMR hierarchy, and the FMG-cycle on a uniform base 
grid. By contrast, the multigrid method presented here 
uses not only V-cyclcs, but also FMG-cycles on the grids 
of the AMR hierarchy. An FMG-cycle on the hierarchical 
grids is also implemented in the approach of Matsumoto & 
Hanawa (2003a), and the same strategy is adopted here. 

Figure 5 shows the grids used in the multigrid cycles 
for the case N x = N y = N z = 8. For convenience, a grid is 
labeled by 

with ee[0,£ max ] and me[0,m max ], (37) 

where £ and m denote the grid- level of the AMR hierarchy 
and the coarsening level, respectively. For example, the 
hatched grid in Figure 5e is labeled as f2 2 , where the grid 
of £ = 1 is coarsened by a factor 2 2 . The coarsening level 
of all the grids shown in Figure 5a is m = 0, while that of 
the grids shown in Figure 56 is m = 1, irrespective of the 
grid-level I. In addition a composite grid is defined as 

Cl m ^fi^U-un^. (38) 

For example, the composite grids shown in Figures 5 a, 56, 
and 5 c are expressed as fi 1 , and f2 2 , respectively. 

An overview of the cycles used by the present numerical 
method is now presented. First, the solution on Ct m for 
< m < m max converges under FMG-cycles, as shown in 
Figures 5 a- 5 c and 6 a. The maximum coarsening level in 
the FMG-cycle is given by m max = log 2 min ( N x , N y , N z ) - 
1. In other words, the grid is coarsened until the number of 
the cells per block is decreased up to two in at least one di- 
rection. At the bottom of the FMG-cycle, which is marked 
by V in Figure 6 a, the solution then converges on VLj 1 for 
< £ < f max and m = m max under V-cycles (the hatched 
grids in Figures 5d-5f). Typically, only one iteration of 



7 

the V-cycle is sufficient, although, for reference, two cy- 
cles of the V-cycle are illustrated in Figure 66. At the 
bottom of the V-cyclc, which is marked by B in Figure 66, 
the solution on the coarsened base grid 17™ for m > m max 
converges under FMG-cycles, as shown in Figures 5<?-5i! 
and 6 c. Finally, at the bottom of the FMG-cycle on the 
base grid, which is marked by E in Figure 6c, an exact 
solution is given according to the boundary conditions, 
since there is only one cell in the computational domain. 

Note that the computation at Cl° of the FMG-cycle ac- 
counts for most of the computational time of the multigrid 
method, because the number of cells in f)° dominates all 
others. The number of cells decreases by a factor of 1/8 
every time the grids are coarsened. This indicates that the 
computation on Q° dominates the overall computational 
cost, and this scheme is scalable to the number of cells 
in a similar way as the hydrodynamics scheme (see also 
Matsumoto & Hanawa 2003a). 

6.2. Smoothing 

The red-black Gauss-Seidel iteration is adopted as a 
smoothing procedure. Conventional red-black Gauss- 
Scidcl iteration can be applied only to a uniform grid. 
This iteration scheme is therefore modified so that it can 
be applied to composite grids, in which grids with different 
resolutions co-exist. 

The discretization of Poisson's equation on the com- 
posite grids is now described. Poisson's equation can be 



expressed as a set of two equations, 

V-g = -^G Pl (39) 

g = -V$. (40) 

These equations are discretized as 

d x 9x,i,j,k + d v g V} i^ tk + d z g z jj tk = -^Gpij,k, (41) 

9x,i+l/2,j,k = -dx^i+l/2,],k, (42) 

9y,i,j+l/2,k — ~dy ^i,j+l/2,fc) (43) 

9z,i,j,k+i/2 = -d z $i,j,k+i/2, (44) 



where g X! i+i/ 2 ,j,k, 9y,i,j+i/2,k, and g z ^j. k +i/2 are the com- 
ponents of gravity defined on the cell surfaces. 

Considering the fine cells adjacent to the interface be- 
tween the fine and coarse cells in Figure 7, the poten- 
tial on the ghost cell, <J> B is required in order to obtain 
9x i-1/2 j k- potential on the ghost cell is given by, 

$ B = bM !±iJ± ] (45) 

$* = 

^ ^ 16— ' (46) 

where $ B is obtained by quadratic interpolation in the 
x-direction, and is obtained by bilinear interpola- 
tion in the y — z plane. The quadratic interpolation in 
the direction normal to the interface satisfies the nec- 
essary conditions for so-called conservative interpolation 
(see Trottenberg, Oosterlee, & Schuller 2001). Using this 
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Fig. 5. Schematic diagram of the coarsening of grids in the multigrid method. The cell-boundaries and block-boundaries are 
denoted by thin and thick lines, respectively, (a-c) Coarsening of grids in the FMG-cycle on the AMR hierarchy. The number of the 
cells per block decreases up to 2 3 . (d—f) V-cycle on the AMR hierarchy. The solution converges sequentially on the hatched blocks. 
(g-i) Coarsening of grids in the FMG-cycle on the uniform base grid. 



procedure, the components of gravity of the fine cells ad- 
jacent to coarse cells are obtained. 

For the coarse cells adjacent to fine cells, the gravity at 
the interface is given by summing the gravity on the cor- 
responding cell surfaces of the fine cells. For the example 
of Figure 7, such a coarse gravity is given by, 

j'+i fc+i 

9x,I+X/2,J,K = 4 Zj5j^,»-l/2,j,fc ( 47 ) 
3=3 k=k 

instead of g^ I+1 / 2 j k = j k- This ^ s a similar 



strategy to the rcfluxing described in § 5.4, and ensures 
that the solution satisfies Gauss's theorem; when the nor- 
mal component of gravity is summed over the surfaces of 
any cell, the sum equals the mass contained by the cell 
multiplied by AnG, 

g-AS = -AirG P AV, (48) 

surface 

where AV denotes the volume of the cell, and AS denotes 
the vector of the cross-section of the cell in the direction 
normal to the cell surface. We also confirmed that rcflux- 
ing the gravity is necessary for second-order accuracy by 
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Fig. 6. Schematic diagram of the multigrid cycles. The lines pointing diagonally downwards from left to right denote the restriction 
operators, while the lines pointing diagonally upwards from left to right denote the prolongation operators. The S symbols denote 
the smoothing operators. The V symbols in panel a denote the V-cycles on the AMR hierarchy shown in panel b. The B symbols in 
panel b denote the FMG-cycles on the uniform base grid shown in panel c. At the point denoted by E, an exact solution is obtained 
according to the boundary conditions. 



a convergence test (see § 8.5). 

This refluxing of the gravity is adopted only for the 
FMG-cyclc on the composite grids shown in Figure 6 a, in 
which the solution converges simultaneously on f2™ over 
the grid-levels t. The refluxing of gravity is not adopted 
in the V-cycle, since the solution converges sequentially 
over the grid- levels, 

According to the discretization of equation (41), the 
Gauss-Scidcl iteration is expressed as, 



(49) 



R 



i,j,k ■ ®x9x,i,j,k $yQy : i,j : k H~ $z9z,i,j,k ~t~ ^^^Pz^i,j,k 

= £'!',.,> + A7rGp zAtjtk , (50) 

where h denotes the cell width, Rij.k denotes the residual, 
and C denotes the Laplacian operator. Equation (49) up- 
dates &i,j,k to ( &"™ at every iteration. Typically only two 
iterations are required for the FMG-cycle, and one itera- 



10 



Tomoaki Matsumoto 



[Vol. 59, 



u...., 



J. 





■ 


1/ 




1 


■ ■ 










■ ■ ■ 


■ ■ 


| 









-J 




■ 






>. 




■ ■ 









Fig. 7. Interface between the fine and coarse cells for the multigrid method. 



tion for the V-cycle. Hereafter, ^ffk = C^s (®i,j,k,Pi,j,k) 
refers to the Gauss-Seidel iteration given by equations (49) 
and (50). 

6.3. Prolongation and Restriction 

The full-weight prolongation is adopted here, and is 
given by 



1 

64 

9f*f 



27$f, 



J,K " 



^?±l,J±l,K±l 



■*I,J,K±l) 



-I±l,J,K 



(51) 



where (I,J,K) indicates a coarse cell overlapping with a 
fine cell {i,j,k), and the sign of ± depends on the parity 
of i, j, and fc. 

The restriction procedure is performed by averaging the 
values of the fine cells (i,j,k) which overlap the corre- 
sponding coarse cell (I,J,K), 



*£^=tf* h = i£££*fe,it- (52) 



i—i j—j k—k 



The prolongation and restriction introduced above is 
used in common with all the multigrid cycles. 



6.. 



FMG- Cycle on AMR hierarchy 



An FMG-cycle based on the standard algorithm of 
the FMG-cycle for linear equations (see e.g., Press & 
Tcukolsky 1991) is implemented on the composite grids 
of the AMR hierarchy, because the basic operators are 
prepared on the composite grids. The smoothing, pro- 
longation, and restriction procedures are given by equa- 



tions (49), (52), and (52), respectively. In the smoothing 
procedure, the refluxing of the gravity is performed ac- 
cording to equation (47). In the prolongation, the vari- 
ables on fl™ are transferred onto Similarly, the re- 
striction procedure transfers variables on fl™ onto 
The overlap of coarse grids on fine grids, f^Hfi^, is 
therefore not taken into account in the FMG-cycle. 

6.5. V -Cycle on AMR hierarchy 

hi the V-cycle, the multilevel adaptive technique 
(MLAT) is adopted by using the full approximation 
scheme (FAS) (see e.g., Trottenbcrg, Oosterlee, & Schuller 
2001). The solution is iterated towards convergence on 
O™, where the boundary condition at <9S1™ is obtained 
from according to equations (45) and (46). The pro- 
longation procedure transfers variables from j to fl™, 
while the restriction procedure transfers variables from 

to n%_ v 

We now describe the method of the V-cycle in terms 
of the fine and coarse grids, SI™ and SVg]_i- As a prc- 
smoothing procedure, the smoothing procedure is applied 
to the initial guess of <J>' 1 on SI™, 

$' 1 



(53) 



The density on Sl" L _ 1 is then obtained as 



p H '' 
where 



If/ 



onn^nn? (54) 

on the remaining part of W[Li 



(55) 



is the so-called r-correction. The initial guess on SlJLi 
given by, 
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H.- 



H 



on Q^nQf 
on the remaining part of Q 



Using p H '* and Q H * , the approximate solution $ H is ob- 
tained from the coarser grids. Then, i> /l is updated to & h 
on nj 1 according to, 
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.(56) 



$' 1 = $' ; 



■if** 



Finally, post-smoothing is applied to $ 

=,/i,ncw 1 



-GS 



(57) 



(58) 



Thus, by this algorithm, the initial guess & h is updated 
to $^ ncw . 

6.6. Utilization of Multigrid Method 

The front-end of the multigrid method described here 
is the FMG-cycle on the AMR hierarchy. Given density 
p, the boundary condition for <!>, and the initial guess <£>°, 
the FMG-cycle is called iteratively, and <fr™ is updated to 
(jjn+i ky means f the following procedure, 



n+1 



^FMG (0) Ri,j,k) ; 



(59) 
(60) 



where £pMG P) denotes the Poisson solver of the FMG- 
cycle on the AMR hierarchy with initial guess $ and the 
right-hand side of Poisson's equation p. The refluxing 
of the gravity is implemented in both the C and £pM G 
operators. The FMG-cycle £p^ G is always called for the 
fixed boundary condition of zero, because the boundary 
condition is transmitted to the residual R according to 
equation (59). 

The solution converges under the iterative cycle of equa- 
tions (59) and (60), reducing the absolute value of the 
residual \R\. In the problem of cloud collapse, the solu- 
tion of $ in the previous time step is adopted as the initial 
guess, and just a single cycle reduces the residual typically 



by \Rh 2 \ m ^/ (4irGph 2 



io- 



7. Parallelization and Vectorization 

The code is written entirely in Fortran90, and paral- 
lelized by the MPI library. The data is partitioned into 
the computational nodes used for the parallel calculation 
in units of the blocks; all cells within a given block are 
assigned to the same node. For a given grid-level, all the 
blocks are ordered by means of the Peano-Hilbert space 
filling curve, and the blocks are then assigned to the nodes 
according to this order. The vectorization in block units is 
computationally intensive. Since a block has N x x N y x N z 
cells, the vector- length is of the order N x x N y x N z . The 
block sizes, N x , N y , and N z are set at the initial config- 
uration, and the choice of the block size depends on the 
machine. In the typical case N x = N y = N z = 8 is chosen, 
and the vector-length is therefore about 512, which is ac- 
ceptable for the vector processor used in the calculations. 
In the next section, convergence tests are performed by 
changing the block sizes, N x , N y , and N z . 



8.1. 

Convergence tests for the linearized MHD waves (fast, 
slow, Alfven, and entropy waves) are shown. Following 
Crockett et al. (2005) and Gardiner & Stone (2005), the 
waves propagation at a slope of 2:1 is performed on uni- 
form grids and AMR hierarchical grids in the x — y plane. 

The unperturbed state is set at, po = 1, Pq = 1, Vq = 0, 

. T 

1 

heat ratio is 7 = 5/3. All the linearized MHD waves have 
wave number of k = 2ir(2, 1,0) T . The wavelength is there- 
fore A = 1/ \/E, and the computational domain includes 
two net waves (Fig. 8). 

The initial perturbations of the fast and slow waves are 
set based on the eigenmode, such as, 



and B = (^751 "TJ'O^ ina;,|/£ [0,1]. The specific 



6p 
<5u|l 
6v± 
SB ± 
SP 

/ 



Po 

CF/S 

\\B±_ c F / s 



iirpo c 



B±- 



(7"1)( 
$p Gr t sin(fc 
1 



Poc F/s 
r), 



F/S 
C F/S 

B ± SB± 



)-( 7 -2)pc2 ) 



-F/S 



2 

4irp a 

iirpo 
Po 
Po 



c 2 A ± 



4c? c' 2 



(61) 
(62) 

(63) 
(64) 
(65) 



where v\\ and v± denote the parallel and perpendicular 
components of v with respect to k. Similarly, B\\ and B± 
denote the parallel and perpendicular components of B 
with respect to k. The subscripts F and S represent the 
fast and slow waves. The amplitude of the perturbation 



is set at <5 



pert 



10 5 . For the Alfven wave, the initial 



perturbation is set as, 



SB Z 



CA 

-B 



J pcrt 



sin(fc 



(66) 



The entropy wave is a simple advection wave, and the 
perturbation is imposed only on the density, 

5p = p 5 pelt sin(fc-r), 

and the velocity of the unperturbed state is set to, 

k 



v 



(67) 



(68) 
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Fig. 8. Initial conditions for propagation of the simple lin- 
earized MHD waves on (a) the uniform grid and (6) the AMR 
hierarchical grid. The rectangles denote the block distribu- 
tions and each block contains N x X N v = 4 2 cells. 



rather than Vo = 0. 

Figure 8 shows the wave patterns at the initial condi- 
tions and the block distribution. The wave propagates in 
the upper right direction. The uniform grid consists of 
4x4 blocks of the base grid £ = (Fig. 8 a). In the AMR 
hierarchical grid (Fig. 86), the region of x < 0.5 is covered 
by fine grids of I = 1. In this test, the AMR hierarchical 
grid is static; the block distribution shown in Figure 8 is 
preserved in the course of the calculations in order to in- 
vestigate wave propagation on the coexistence of the fine 
and coarse cells. For both the uniform and AMR grids, 
the periodic boundary condition is imposed on the x = 0, 1 
and y = 0, 1. The three-dimensional code is applied to this 
problem even though it has a two-dimensional symmetry, 




Fig. 9. Li norm of error as a function of cell width of 
the coarsest grid £ = 0, ftbasei f° r ( a ) the left side (x < 0.5) 
and (6) the right side (x > 0.5) of the computational domain. 
The lines with filled circles, open circles, filled diamonds, and 
open diamonds denote errors for slow, fast, Alfven and en- 
tropy waves, respectively. The solid and dashed lines are for 
the AMR and uniform grids, respectively. The dotted lines 
indicate the relationship of errors in proportion to , and 

that in proportion to /ibasc- 

with all the variables being constant in the z-direction. 

We performed the convergence test by changing the cell 
number inside a block as 4 < N x ,N y < 256 for the uniform 
grids, and 4 < N x , N y < 128 for the AMR grids. This 
corresponds to a change in the resolution from h = 1/16 
to 1/1024 for the uniform grids, from h = 1/16 to 1/512 
for grid-level £ = of the AMR grids, and from h = 1/32 to 
1/1024 for grid-level £ — I, respectively, where h = Ax = 
Ay. Each run is halted when a wave is propagated by a 
distance of three wavelengths so that all the waves sweep 
the entire computational domain. A Li norm of the error 
is estimated by 



E 



u 



u 



m,0 



AS, 



(69) 
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Fig. 10. Decay rate as a function of cell width for slow, 
fast, Alfven, and entropy waves in the uniform grid. The 
lines with filled circles, open circles, filled diamonds, and open 
diamonds denote errors for slow, fast, Alfven and entropy 
waves, respectively. The dashed lines indicate the relationship 
of errors in proportion to /i 3 , and that in proportion to h 2 . 

where {t/™} = Uij denotes the conservative variables de- 
fined by equation (12) at the last stage, |^™'°} = 
denotes that of the initial state, and AS^.j = Ax% jAy% j 
denotes the cell surface. For the case of the AMR grids, 
only the cells of I = 1 are taken into account for the sum 
in the region of x < 0.5 since grid-levels of i = and 1 
overlap in this region. 

Figure 9 shows the L\ norm for regions of x < 0.5 and 
x > 0.5, separately. The L\ norms of the uniform grid 
(dashed lines) in Figures 9a and 9b coincide with each 
other because of the symmetry. All the waves exhibit 
second-order accuracy for /ibase ^ 1/64, and first-order ac- 
curacy for ft-baso ^ 1/32, irrespective of the grid types and 
the wave modes, where /ibase denotes the cell width of the 
coarsest grid (the base grid). This indicates that the nu- 
merical method basically has second-order accuracy, but 
the grids of /ibase ^ 1/32 are too coarse to resolve the 
waves. 

The AMR girds always exhibit smaller L\ norms than 
the uniform grids for given a wave mode and /ibase- This 
is because the AMR grids resolve the waves using a grid in 
the region of x < 0.5 that is two times finer than the cor- 
responding uniform grids. This indicates that the AMR 
improves the solution of the wave propagation for all the 
wave. However, it should also be noted that the solution 
is improved more when the resolution of the uniform grid 
is increased by a factor two. Comparing the left and right 
sides for the AMR grids (Figs. 9a and 96), the left region 
shows slightly smaller errors than the right region for the 
fast wave because the left side has finer resolution. In 
contrast, the left side has slightly larger errors than the 
right side for the slow and entropy waves, in spite of the 
finer resolution of the left side, indicating the effects of 
reflection of the waves at the boundaries between the fine 
and coarse grids. For the Alfven wave, the significant dif- 



ference is not observed between the left and right sides. 

For the AMR girds, a wave propagates from the fine 
grid into the coarse grid at x = 0.5, and it also propagates 
from the coarse grid to the fine grid at x = and 1 via 
the periodic boundary condition. A significant reflection 
of the waves at x = 0.5 is observed in the cases of very 
coarse grids (/ibase ^ 1/32) in the v y component of the fast 
wave and the v x component of the slow wave. However, 
it should be noted that the amplitude of the components 
showing reflection is much smaller than that of the other 
components, e.g., v y has an amplitude of only 10 % of v x 
for the fast wave. 

Estimation of decay rate of waves is important for il- 
lustrating properties of a numerical method as indicated 
by Crockett et al. (2005), and we also estimate the decay 
rate for the test of the wave propagation on the uniform 
grids described above. The amplitude of the wave of the 
iTi-th mode is estimated by 

w m =Y,\l m (U i , j -U )\ (70) 

where Uij denotes the conservative variable, Uq denotes 
that in the unperturbed state, and l m denotes the left 
eigenvector of a wave mode evaluated in the unperturbed 
state. The superscript m indicates the modes of the MHD 
waves: the fast, slow, Alfven, and entropy waves. 

Figure 10 shows the decay rates for all the MHD waves 
as a function of the cell width. The decay rates for all the 
wave is almost in proportion to h 3 for h < 1/ 64 and h 2 for 
h > 1/64. This dependence of the decay rate on the cell 
width is consistent with Crockett et al. (2005), indicating 
that the region of the second-order accuracy exhibits the 
third-order accuracy in the decay rate. The third-order 
decay rate of the scheme presented here is also reported by 
Sugimoto et al. (2004) (see their Fig. 1), in which the same 
numerical flux presented here is adopted. In contrast, 
Ryu et al. (1995) reported a second-order decay rate by 
using a TVD scheme. The TVD approach is also adopted 
here but is based on the scheme of Fukuda & Hanawa 
(1999), where the MUSCL extrapolation is applied to the 
amplitude of the each eigenmode l m U rather than the 
primitive variables Q. 

8.2. Magnetic Flux Tube 

The magnetic flux tube problem is proposed by 
Crockett et al. (2005), and is a test for stiffness of a MHD 
scheme. The model set-up is the same as that of Crockett 
et al. (2005). The unperturbed state is set as p = 1, v = 0, 
B x = 0, and B z = in x,y € [—0.5,0.5]. The specific heat 
ratio is 7 = 5/3. The computational domain is divided into 
a magnetized region of x G [—0.2,0.2] and non-magnetized 
regions. The magnetized region and the non-magnetized 
regions are balanced via the thermal pressure and mag- 
netic pressure; B y = V807T and P = 1 in x € [—0.2,0.2], 
and By = and P = 11 in the remaining regions. The 
periodic boundary condition is imposed on the x = ±0.5 
and y = ±0.5. On this unperturbed state, the sinusoidal 
transverse velocity of v x = <5 pe rtCASm(27ry) is imposed, and 
the amplitude of the velocity is set at <5 per t = 0.01. This 
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Fig. 11. Magnetic flux tube problem calculated by various models for (top left) kinetic energy, (top right) internal thermal energy, 
(bottom left) magnetic energy, and (bottom right) L\ norm of V ■ B. The abscissa denotes the time normalized by Alfven crossing 
time t a = l/c A = 0.2236. 



Table 1. Parameters for Magnetic Flux Tube 



model 


orientation 


cell width (h) 


grid 




AU512 


aligned 


1/512 


uniform 




AAMR 


aligned 


1/256- 1/1024 


AMR (£ m ax 


= 2) 


IU256 


inclined 


a/2/256 


uniform 




IU512 


inclined 


a/2/512 


uniform 




IU1024 


inclined 


^2/1024 


uniform 




IAMR 


inclined 


V2/256- V2/1024 


AMR (£ max 


= 2) 



model is calculated using two types of grid: a uniform grid 
of h = 1/512 (model AU512), and an AMR grid (model 
AAMR), as shown in Table 1. In the AMR grid, the 
maximum grid- level is set as £ max = 2, and the following 
refinement criterion is adopted, 

max[£ {pi,3,k),£(Pi,3,k)\ > KT 1 , (71) 

£{Qi,j,k) = 1 • (72) 

9i,j,k 



This criterion captures the curvature of the density and 
pressure profiles, and the refinement is attributed mainly 
to the pressure jumps in this problem. Models with uni- 
form grids are calculated until t = 6 (= 26.83t a ), and the 
models with the AMR grid are halted at t = 0.2236 (= t a ) 
because of computational cost, where t a = 1/ca = l/y/20 
denotes Alfven crossing time. 

Grid-inclined flux tube models are also examined. The 
model set-up is the same as that of the grid-aligned model 
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Fig. 12. Magnetic flux tube for models (a) AAMR and (6) IAMR. at the stage of t/t a = 1. The gray scales denote the velocity 
|u|, and the white rectangles denote the block distributions, (a) The region of x £ [—0.2,0.2] is magnetized, and (b) the diagonal 
region is non-magnetized. The boundaries between the magnetic and non-magnetic regions are resolved by the fine blocks. 




Fig. 13. Magnetic flux tube problem calculated by the uniform grid models for (top left) kinetic energy, (top right) internal energy, 
(bottom left) magnetic energy, and (bottom right) L\ norm of V ■ B. The color legend is same as that of Figure 11. 
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described above, but the coordinates are rotated at an 
angle of 45°. Because of this rotation, the computational 
box is larger than the grid-aligned model by a factor s/2: 
x,y £ [I/a/2, — l/v^]- We introduce a taper region at the 
boundaries between the magnetic and non-magnetic re- 
gions, in order to reduce the initial V • B error. At the 
taper region, the pressure and magnetic field have a tanh 
profile, the scale length of which is twice the cell width. 
For the inclined cases, four models are constructed as 
shown in Table 1; three models arc calculated using a 
uniform grid by changing the resolution from h = -\/2/256 
to \/2/1024 (models IU256, IU512, and IU1024), and one 
using the AMR grid (model IAMR). 

Figure 11 shows the kinetic energy, internal thermal 
energy, magnetic energy, and L\ norm of V • B error as 
a function of time for < t < t a . The three former en- 
ergies are normalized using the initial values in the plots. 
The divergence error is calculated as follows (see the ninth 
component of equations [12]-[14]): 



(V-B)^ /2 



(73) 



In the period < t < t a , all the models are calculated 
without any indications of the instability. The two grid- 
aligned models (AU512 and AAMR) exhibit good agree- 
ment with each other, and L\ norms of the V • B are 
remain low. The block distribution of the AMR mod- 
els is shown in Figure 12 a, and the fine blocks cover the 
boundaries between the magnetic and non-magnetic re- 
gions. For the grid-inclined models, the three models 
(IU256, IU512, and IU1024) are calculated with uniform 
grids of different resolutions, and these models show ap- 
proximately first-order convergence (cx h 18 ). The model 
with the AMR grid (IAMR) shows an evolution very sim- 
ilar to the model with the finest uniform grid (IU1024), 
as shown in Figrue 11. For model IAMR, the boundaries 
between the magnetic and non-magnetic regions are re- 
solved by the finest blocks whose cell width is the same as 
that of model IU1024 (Fig. 126). This implies that a domi- 
nant error is caused in the boundary between the magnetic 
and non-magnetic regions. However, model IAMR has a 
V • B error as large as that of model IU246, which has the 
same resolution as the coarsest grid-level of model IAMR. 
Examining the velocity distributions shown in Figure 12, 
model AAMR exhibits good agreement with IAMR when 
it rotates at angle of 45° ; note that the diagonal region of 
Figure 126 corresponds to regions of x < —0.2 and x > 0.2 
in Figure 12 a. 

Figure 13 is the same as Figure 11 except for < t < 6 
and the uniform grid models. The long time calculations 
are performed successfully without any instabilities for 
all the uniform gird models. The kinetic energy of the 
grid-aligned model (AU512) decays slightly. For the grid- 
inclined models, the model with higher resolution shows 
less decay in the kinetic energy, but the kinetic energies of 
the finest grid-inclined model (IU1024) is about half those 
of the grid-aligned model (AU512) at the peak of the final 
osculation. The internal energies show different tenden- 



cies in the grid-aligned and -inclined models, but their 
gradual increases are within the order of 0.1%. The L\ 
norms of the V • B remain low in the grid-aligned model 
(model AU512) and the grid-aligned model with the fine 
grid (model IU1024). 

8.3. Double Mach Reflection 

We consider the double Mach reflection problem as a 
hydrodynamics test problem. This test problem was ini- 
tially proposed by Woodward & Colella (1984), and has 
since been widely used for testing high resolution schemes. 
A planar shock of Mach 10 travels in a medium of p = 1.4, 
P = 1, and 7 = 1.4 with incident angle of 60° against a 
rigid wall. The incident shock interacts with the wall, 
and a complicated structure develops featuring a strong 
and weak reflected shock, contact discontinuities, and a 
small jet at the wall (lower boundary; 1/6 < x < 4 and 
y = 0). The computational domain given by x € [0,4] and 
y £ [0, 1] is covered by 32 x 8 blocks at £ = 0, and the 
maximum grid-level is set at £ max = 4. Each grid has 8 2 
cells (N x = N y = 8), and so the minimum and maximum 
resolutions are h = 1/64 and 1/1024, respectively. The 
following refinement criterion is adopted, 



max[£ (pi,j,k),£(Pij,k)] > 10" 



(74) 



where £ is defined by equation (72). Three-dimensional 
code is applied to this problem even though it has a two- 
dimensional symmetry, with all the variables kept con- 
stant in the z-direction. In order to maintain this two- 
dimensional symmetry, the condition F z = is imposed. 
The Roe scheme is modified here based on Kim et al. 
(2003) in order to avoid the carbuncle instability at shock 
waves. The carbuncle instability refers to a protuber- 
ant shock profile, which is an unphysical phenomenon and 
mashes up shock waves (see, e.g., Perry & Imlay 1988). 
This is thought to be attributed to a low numerical dif- 
fusion of the Roe's scheme, and calculations having high 
resolutions tend to suffer from this instability. 

Figure 14a shows the density distribution at t = 0.2. 
The shock waves are resolved by finest grids as shown 
in Figure 146. Figure 14c shows an enlargement of 
Figure 14a around the double Mach stems. An eddy 
is clearly resolved at the end of the jet, and it can be 
seen even more clearly in the entropy distribution Pj p 1 
(Fig. 14d). The shock wave is resolved clearly at x = 2.8 
because the carbuncle instability was eliminated. Without 
eliminating the carbuncle instability, a serious instability 
was observed in this shock wave, and the eddy was also 
mashed up. At (x,y) — (2.3,0.1), the discontinuity of the 
sheer begins to twist, and another eddy may develop as 
shown by Shi et al. (2003). 

The double Mach reflection problem is also calculated 
using uniform grids changing the resolution from h = 1 /64 
to 1/1024 in order to confirm the convergence. This test 
exhibits an order of convergence of 1.23, which is slightly 
higher than the first order, in spite of implementing the 
second order accuracy. This may be because the slope 
limiter in the MUSCL extrapolations reduces the accuracy 
to the first order in the regions having shock waves and 
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Fig. 14. Double Mach reflection problem at t = 0.2, (a) density distribution in [0,3] X [0,1], (b) block distribution in [0,3] X [0,1], 
(c) enlargement of Figure 14a ([2,2.9] X [0,0.5]), and (d) distribution of P/ p 1 ■ In (a) and (c) 30 contour lines are shown in the range 
1.5 < p < 22.9705, while in (d) 22 contour lines are shown in the range 1.0 < P/ p 1 < 11.5. 



contact discontinuities. 

84. MED Rotor Problem 

The MHD rotor problem was first proposed by Balsara 
& Spicer (1999), and it tests the propagation of non- linear 
Alfven waves. The model set-up is in the same manner as 
that of Toth (2000) and Crockett et al. (2005). The com- 
putational domain is a unit square of x,y £ [0,1]. In the 
initial stage, a uniform cylinder with p = 10, P = 1, and 
radius 0.1 rotates at an angular velocity of 20. The ambi- 
ent medium is at rest with p=l, P = l, and v = 0. At the 
boundary between the cylinder and ambient medium, a ta- 
per region is used in order to reduce the initial transition 
(see Toth 2000). The computational domain is subject 
to a uniform magnetic field, (B x ,B y ,B z ) = (5,0,0). The 
adiabatic index of the gas is 7 = 1.4. The computational 
domain is covered by 16 x 16 blocks at £ = 0, and each 
block consists of 8 x 8 cells (N x = N y = 8). The maximum 
grid-level is set at f max = 2. The coarsest and finest grids 
therefore exhibit an effective resolution of h = 1/128 and 
1/512, respectively. The following refinement criterion is 
adopted, 

rnax^p^, k ),£(Pu, fc )] > 1, (75) 

£(Qu,k) = [£l{qi,j,k] + £ 2 y {<ii,i,k) + £ 2 A% j ,k)} 1/2 , (76) 

h 2 dlqt,i,k ^ 

h\d x q l+1/2 ,j,k\ +fi>\d x q i - 1 /2 ! j,k \ + e^S^ij.fc ' 
and £ y and £ z are defined in a similar manner to £ x , where 
e = 10~ 2 (e.g., Fryxell et al. 2000). 

Figures Iba-lbd show the density, thermal pressure, 
Mach number, magnetic pressure at t = 0.15. The finest 



grids capture the outer shock fronts and inner complex 
structures, as shown in Figure 15/. The distribution of the 
physical variables plotted here exhibit excellent agreement 
with Figure 18 of Toth (2000) and Figure 12 of Crockett 
et al. (2005). Figure 15 e shows the amplitude of the mag- 
netic divergence error normalized by the cell width and the 
local magnetic field strength, h(V ■ B)/\B\, where V • B 
is estimated using equation (73). The divergence error 
reaches a maximum value of 1.5 x 10 -2 at the inner dis- 
continuity, and is 2 x 10 -3 at the outer shock fronts. The 
divergence cleaning of Dedner et al. (2002) keeps the di- 
vergence error small. 

8.5. Accuracy of Multigrid Method 

The accuracy of the multigrid method is examined us- 
ing the approach of Matsumoto & Hanawa (2003a). Two 
uniform spheres of masses 1 and 2, both having a ra- 
dius of 6/1024, are located at (x,y,z) = (12/1024,0,0) 
and ( — 12/1024,0,0) in the computational domain x,y,z £ 
[—0.5,0.5]. A Dirichlet physical boundary condition is 
imposed on <!>, and it is evaluated by the multi-pole ex- 
pansion of the density distribution at £ = 0, where the 
monopole, dipolc, and quadruple moments are taken into 
account. 

The computation was started from the initial guess $ = 
0. The residual decreased by more than a factor of several 
hundred with each FMG iteration given by equations (59) 
and (60). After 10 iterations, the residual was of the order 
of the round-off error. The numerical solution for the 
gravity obtained by this computation is compared to the 
exact analytic solution. 

Figure 16 shows the relative error of the numerically 
calculated gravity, \g — g ex | / |«7 ex |, on a logarithmic scale, 
where g cx denotes the exact gravity obtained analytically. 
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Fig. 15. MHD rotor problem solved by the AMR code, (a) Density, (6) thermal pressure, (c) Mach number, (d) magnetic pressure, 
(e) magnetic monopole, and (J) grid distribution are shown at t = 0.15. In (a)-(d) 30 contour lines are shown in the ranges of 
0.483 < p < 12.95, 0.0202 < P < 2.008, < \v\/c a < 8.18, and 0.0177 < B 2 /8tt < 2.642, respectively. 
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Fig. 16. Distribution of the relative error of the numerically 
calculated gravity in the z = plane with N x = N y = N z = 16 
after 10 iterations of the FMG-cycles. The gray scale denotes 
log ( \g — g cx | / | g ex | ) , and the rectangles indicate the block dis- 
tribution. Panel b is an enlargement of panel a. 

In this calculation, N x = N y = N z = 16 and £ max = 4 
are adopted. The block distribution is also shown in 
Figure 16. The maximum error occurs at the edges of the 
spheres due to the discretization error that arises from 
the sharp density contrast of the spheres. In the other do- 
mains, the error is less than ~ 10 -3 . It is noteworthy that 
no significant error appears at the interfaces between the 
coarse and fine grids. This is attributed to the refluxing 
of the gravity described in § 6.2. 

Figure 17 shows the dependence of the error on resolu- 
tion. The error of the gravity is measured by changing the 
number of cells inside a block, N x , N y , and N z , but fix- 
ing the block distribution. The ordinates of Figures 17a, 
176, and 17c denote the errors measured by the average, 
root mean square, and maximum values, corresponding 
to the L\, L2, and norms, respectively. The abscissa 
denotes the cell width of the coarsest grid I = 0. The 
five lines show the errors at the grid-levels of I = 0, •■■,4, 



while the five points in the lines denote the errors for 
N x =N y = N z = 2, 4, 8, 16, and 32. All the lines, except 
for £ = 4 in Figure 17 c, exhibit second-order accuracy for 
the multigrid method presented here. By contrast, the 
line of I = 4 in Figure 17c shows a first-order accuracy. 
This is attributed to the discretization error of the den- 
sity imposed on this grid-level. 

8. 6. Convergence of Multigrid Method 

The reduction of the residual defined by equation (50), 
is measured to evaluate the efficiency of the iteration for 
the same problem as described in § 8.5. Figure 18 shows 
the maximum residual in each grid-level as a function of 
the number of the FMG-cycles defined by equations (59) 
and (60) for the cases N X = N V = N Z = 8, 16, and 32. The 
residuals plotted here are multiplied by h 2 so that they 
have a dimension of \I>. For the case N x = N y = N z = 8 
(Fig. 18a), the residuals in all the grid- levels decrease in 
proportion to exp(— 6n). After 7 iterations of the FMG- 
cycles, the residuals at all grid-levels reach the round-off 
error. On the other hand, for the cases N x = N y = N z = 16 
and 32 (Figs. 18& and 18c), the rate at which the residuals 
decrease on the coarse grid-levels are slower than those of 
the fine grid-levels. The slower convergence is due to the 
boundary condition imposed on the coarsest grid. 

8. 7. Collapse and Fragmentation of an Isothermal Cloud 

The present numerical technique is applied to the prob- 
lem of the collapse and fragmentation of an isothermal 
cloud as a gravitational hydrodynamics test problem. 
While isothermal collapse has been calculated by many 
authors, the particular model given in Bate & Burkcrt 
(1997) and Truelovc et al. (1998) is followed here. The 
initial cloud has a uniform, spherical density distribution, 
and rotates rigidly around the z-axis. The mass of the 
cloud is IMq, and the radius is i? c = 5 x 10 16 cm, and so the 
density of the cloud is therefore po = 3.79 x 10~ 18 gcm~ 3 . 
Conventionally, a cloud is characterized by two global 
quantities, a = E th /\E grav \ and (3 = £ , r ot/|-E , grav| 1 where 
-Ethi -Erot, and E grav denote the thermal, rotation, and 
gravitational energies. The cloud here has a = 0.26 and 
f3 = 0.16. The sound speed and the angular velocity are 
obtained as c s = 0.166kms _1 and £1 = 7.14 x 10~ 13 s -1 , 
respectively. The cloud is perturbed by a bar perturba- 
tion with an amplitude of 10%; p = po(l + 0.1cos2</>). The 
cloud is embedded in an ambient gas, whose density is 
O.Olpo- 

The computational domain is x,y,z G [— 2R C , 2R C , ] x 
[—2R C ,2R C ,] x [0,2i? c ,], and a mirror boundary condition 
is imposed on the z = plane. Blocks of N. x = N y =N Z = 8 
are adopted. At the initial stage, 16 x 16 x 8 blocks are 
distributed at the grid- level of £ = 0. The cloud radius R c 
is therefore resolved by 32 cells. This initial resolution is 
same as that of Truelove et al. (1998). The Jeans condi- 
tion is employed as a refinement criterion; the blocks are 
refined when the Jeans length is shorter than 8 times of 
cell width; (ncl/Gp) 1 ^ 2 < Sh. This refinement criterion is 
twice as severe as that of Truelove et al. (1997). 

Figure 19 shows the maximum density as a function of 
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Fig. 17. The relative error of the numerically computed 
gravity \g — g ex | / |g ox | as a function of /ibasc (cell width of 
the coarsest grids). The error is measured by (a) average, 
(6) root mean square, and (c) maximum on each grid-level 
separately. The open diamonds, open circles, filled diamonds, 
filled circles, and filled square denote the errors on grids of 
1 = 0,1,2,3, and 4, respectively. The dashed lines indicate 
the relationship of (a-c) errors in proportion to , and (c) 
those in proportion to /ibase 



Fig. 18. The maximum residual (|/i 2 -R|max) as a func- 
tion of the number of the FMG-cycles for the cases (a) 



8, (b) N x 



N v 



16, and (c) 



N x = Ny = N z 
N x = N y = N% = 32. The open diamonds, open cir- 
cles, filled diamonds, filled circles, and filled square de- 
note the residuals measured on grids of I = 0,1,2, 3, and 
4, respectively. The dashed line displays the relationships 
|/l 2 -R n | ocexp(— An) and exp(— 6n), where n denotes the 
number of the FMG-cycles. 
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Fig. 19. Maximum density as a function of time for the 
problem of the collapse and fragmentation of an isothermal 
cloud. 

time. The central density initially increases gradually as 
the cloud collapses. A thin disk forms due to the fast 
rotation of the cloud when the central density reaches a 



value of around p max = 2 x 10 15 gcm 3 at t — 1.25 x 10 12 s 
(= 1.16%), where t s = (3tt /ZZGpo) 1 / 2 = 5.60 x 10 12 s, and 
this disk bounces in the z-direction. After the bounce, 
the density increases more rapidly, and approaches a sin- 
gularity at t = 1.35 x 10 12 s (= 1.25 

Figure 20 a shows the density distribution in the z = 
plane when p max = 1.02 x 10 5 po- The disk has two density 
peaks, which are resolved by the fine blocks. Each density 
peak has a filamentary structure. This density structure 
was also found by Truelove et al. (1998) (see their Fig. 12). 

Figure 206 shows the upper left bar whenp max = 1.01 x 
10 8 po- The bar collapses to form a very narrow filament. 
This is also consistent with Truelove et al. (1998) (sec their 
Fig. 13). 

8.8. Collapse and Outflow Formation of a Cloud Core 
with Slow Rotation and Oblique Magnetic Field 

The collapse of a cloud with magnetic field parallel to 
the rotation axis has been simulated by Machida et al. 
(2004); Machida et al. (2005a); Machida et al. (2005b); 
Zicglcr (2005); Banerjee & Pudritz (2006); Fromang et 
al. (2006), while the collapse of an oblique magnetic field 
has been simulated by Matsumoto & Tomisaka (2004); 
Machida et al. (2006). In the case of an oblique magnetic 
field, the rotation axis changes its direction due to the 
anisotropic magnetic braking during the collapse, and an 
outflow is ejected after an adiabatic core formation. 

In this paper, the model MF45 of Matsumoto & 
Tomisaka (2004) was calculated using the AMR presented 
here. The initial cloud has the density profile of a critical 
Bonnor-Ebert sphere (Ebert 1955; Bonnor 1956), but the 
density is increased by a factor of 1.68. The central den- 
sity is po = 1 x 10~ 19 gcm~ 3 , and the radius of the cloud 
is R c = 5.49 x 10 17 cm. The cloud rotates slowly at an an- 
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Fig. 20. Collapse and fragmentation of a rotating isother- 
mal cloud. The gray scale denotes logarithmic density in the 
z = plane, and the rectangles indicate the block distribution, 
(a) The central region of the cloud when p max = 1-02 X 10 5 po 
on the grid-levels i = 3 — 7. (b) Enlargement around the up- 
per left fragment when p m ax = 1.01 X 10 8 po on grid-levels 
£ = 6-12. 

gular velocity of Qq = 7.11 x 10 _7 yr _1 (= O.lbt^ 1 ), where 
the freefall time is ts = 2.10 x 10 5 yr. The initial mag- 
netic field is inclined at an angle of 45° with respective 
to the z-axis (rotation axis), and has a strength 18.6 pG. 
The cloud has parameters given by a = 0.5, (3 = 0.02, and 



J\E e 



■■ 0.721, where E mag denotes the total mag- 



netic energy inside the cloud. 

In order to mimic the formation of the first stellar core, 
the equation of state is changed according to the density: 
an isothermal gas of T = 10 K (c s = 0.19kms _1 ) is assumed 
when density is less than p CI = 2 x 10 _13 gcm~ 3 , while a 
polytrope gas of 7 = 5/3 is assumed when density is higher 
than p cr . 

Figure 21a shows the central region of (287 AU) 3 at 
t = 4.81 x 10 5 yr. As calculated by Matsumoto & Tomisaka 
(2004), the outflow is ejected from the region near the 
first core. The outflow speed reaches 10c s , which is also 
consistent with a previous calculation. 
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Fig. 21. (a) Three-dimensional view of protostellar collapse 
and outflow formation at t = 4.81 X 10 s yr. The 6 disk-shaped 
isosurfaces are shown for 4.7 < logp/ po < 7.0, and the 4 bi-po- 
lar isosurfaces are shown for 2.0 < v r /c a < 9.0. The tubes in- 
dicate the magnetic field lines. The coordinates are shown in 
the units of e s /(47rGpo) 1/ ' 2 - (b) Same as panel abut the block 
distribution is overplayed. The grid-levels is ranged from 11 
to 13. 

Figure 21 & shows the block distribution, where grid- 
levels of i = 11 — 13 are shown. Each grid- level has 8 3 
blocks. This grid distribution reproduces effectively the 
same cell distribution as a nested grid including 64 3 cells 
in each grid-level. 

9. Summary 

SFUMATO, a self-gravitational MHD code applying 
the AMR technique is presented. The grid is configured 



in a block structure. 

The MHD scheme is implemented so that it is of second- 
order spatial accuracy by means of the TVD approach. 
The upwind numerical flux is obtained by the linearized 
Riemann solver. The scheme is fully cell-centered, and 
the divergence error of the magnetic fields is cleaned. The 
convergence tests of the linearized MHD waves exhibit a 
second-order accuracy, and the decay rates shows a third- 
order accuracy. Stiffness of the scheme is confirmed by 
the MHD flux tube problems. 

The self-gravity is solved by a multigrid method com- 
posed of: (1) a FMG-cyclc on the AMR hierarchical grids, 
(2) a V-cycle on these grids, and (3) a FMG-cycle on 
the base grid. The FMG-cycle on the AMR hierarchical 
grids enables a scalable dependence on the number of cells. 
The multigrid method ensures that the solution converges 
rapidly; the residual is reduced by a factor of 10~ 3 — 10~ 2 
every iteration. The multigrid method exhibits a second- 
order spatial accuracy. Moreover, no spurious features ap- 
pear at the interfaces between fine and coarse grid- levels, 
due to the flux conservation at the interface. 

The MHD scheme and multigrid method are combined 
so as to have second-order temporal accuracy. The time- 
marching has two modes: synchronous and adaptive time- 
step modes. The former mode is adopted for problems 
including self-gravity, while the latter mode is used for all 
other problems. 

The AMR code was tested by considering test prob- 
lems given by double Mach reflection, the MHD rotor, 
fragmentation of an isothermal cloud, and outflow forma- 
tion in a collapsing cloud. The present results are in good 
agreement with those of previous studies conducted by the 
present author as well as other authors. 
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